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We study the spin and charge fluctuations of the extended Hubbard model (EHM) with on- 
site interaction U and first neighbor interaction V on the two-dimensional square lattice in the 
weak to intermediate coupling regime. We propose an extension of the two-particle self-consistent 
(ETPSC) approximation that includes the effect of functional derivatives of the pair correlation 
functions on irreducible spin and charge vertices. These functional derivatives were ignored in 
our previous work. We evaluate them assuming particle-hole symmetry. The resulting theory 
satisfies conservations laws and the Mermin- Wagner theorem. Our current results are in much 
better agreement with benchmark Quantum Monte-Carlo (QMC) results. As a function of U and 
V, we can determine the crossover temperatures towards renormalized classical regimes where either 
spin or charge fluctuations dominate. The dominant wave vector is self-determined by the approach. 

PACS numbers: 71.10.Fd, 05.30.Fk, 71.10.-w 

& ' I. INTRODUCTION 

fa' 

, The single-band Hubbard model, with on-site interaction U and nearest-neighbor hopping t, is the simplest model 
of correlated electrons. It has been extensively studied, for example, in the context of high-temperature superconduc- 
tors. Charge ordering phenomena in strongly correlated systems are also quite common, for example in manganites, 
vanadates, cobaltates and various organic conductors.—^ In these cases, and even in high temperature superconductors 
where screening is not perfect, it is useful to investigate the Extended Hubbard model (EHM) that includes the effect 
r"| i of nearest-neighbor Coulomb repulsion V. In that model, charge and spin fluctuations compete over much of the 
O ' phase diagram. 

O , Even at weak to intermediate coupling, (U and V less than the bandwidth W — 8t in d = 2) very few many-body 
methods can treat the Hubbard model while adequately taking into account the Mermin- Wagner theorem. The 
H Fluctuation-Exchange Approximation (FLEX) 3 is such a method. The Two-Particle Self- Consistent approach 
£> ■ (TPSC)^£&£ is an alternate non-perturbative approach that compares much better with benchmark Quantum Monte 
Carlo (QMC) simulations. In addition, contrary to FLEX, it reproduces the pseudogap phenomenon observed in 

r-j : QMCSA 

In this paper, we present an extension of the TPSC approach that allows one to treat the EHM. The accuracy 
of the new approach is checked by comparing with QMC simulations^ The competition between charge and spin 
' fluctuations in parameter space is also explored for two fillings. Whenever possible we compare with results of other 
approaches available in the literature. Studies for the specific parameters appropriate for given compounds will be 
""j— ? ■ the subject of separate papers. The purpose of the present one is to give a detailed description of the theoretical 
approach, which is an extension of the TPSC approach. 

Briefly speaking, the TPSC method was derived from a functional approach by approximating the vertex functions 
as local functions of the imaginary time and space (in other words irreducible vertices independent of frequency and 
wave vector). This leads to irreducible spin and charge vertices that contain two unknown constants, namely the 
on-site pair correlation function and its functional derivative, which can be self-consistently determined from two 
£j \ sum-rules on the structure functions. Unfortunately the number of unknowns becomes more than the number of 
available sum-rules in the EHM, so one has to use another approach. In our previous workjii we simply neglected 
the contribution of functional derivatives of the pair correlation functions. We obtained good agreement with QMC 
results for small values of V when spin fluctuations are still dominant, but this approach did not work so well in the 
opposite limit when charge fluctuations becomes dominant. This is mainly due to ignoring the functional derivative 
of the pair correlation functions. In this paper, we demonstrate the importance of these functional derivative terms, 
explain the way they enter in the formalism and propose a way to use particle-hole symmetry to evaluate them. 

The EHM has been studied within different approaches and approximations. We can only give a sample of the 
literature. Many papers focus on the quarter-filled case on the square lattice or on ladders2 i 12 ' 13 i 14 i 15 i 16 i 17 i 18 . The 
competition between charge and spin orders has also been studied in other dimension: one 19 i 20 i 21 i three— and 
higher dimensions^ 3 , at various fillings. In two dimensions, there are several studies at strong couplin g 24 i 25 ' 26 or on 
lattices different from the square latticej 27 ! 28 The combined effect of charge fluctuations in addition to the usual spin 
fluctuation in favouring one type or another of unconventional superconductivity has also been studied using this 
model 2 ^ 3 ^ 3 ^. 

On the two-dimensional square lattice, QMC studies 10 will be particularly relevant as mentioned above. We have 
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not found any systematic comparison between FLEX approximation for the EHMS 9 . and QMC. The mean field (MF) 
approximation 3 ^ and the renormalization group (RG) approach 3 -^, as well as exact diagonalization^, perturbation 
theory about ordered states 3 - 7 - and trial correlated wave functions^ have been applied to this problem. One of the 
results that we will find at half-filling is a crossover line at U « 4V from dominant antiferromagnetic fluctuations to 
dominant checkerboard charge fluctuations. This is in agreement with several studies^ ' 10 ' 29 ' 36 ' 37 ' 38 as we will discuss. 
Our results hold at the crossover temperature. Most other approaches focus on the ground state. The result U — zV, 
where z is the number of nearest-neighbors, seems to hold quite generally in various dimensions ^ 19 ' 39 Some MF 
results 3 -^ give the possibility of a coexistence region between charge and spin density waves. The Pauli principle 
prevents charge and spin fluctuations to increase at the same time in the paramagnetic regime at finite temperature. 
This is satisfied within our approach. 

In the following we present the formalism, which first focuses on finding a set of closed equations for the irreducible 
vertices. We explain how the functional derivative of the pair correlation functions enter in the evaluation of irreducible 
vertices. In the second part of the theory section, we explain how particle-hole symmetry can be used to estimate the 
functional derivatives and we give the details of the new resulting Extended TPSC (ETPSC) approach. In Sec. IIII| 
we present numerical results and compare them with available QMC results to gauge the accuracy of the approach. 
The last part of that section presents a systematic study on the square lattice of the crossover to the renormalized 
classical regimes for charge and spin fluctuations for two densities, n = 1 and n — 0.75. Finally, in the concluding 
section we discuss the possibility to use the same type of approximation for other systems and point towards possible 
future developments and improvements. 



II. THEORY 



A. Some exact results, factorization and approximation for the functional derivatives 

The extended Hubbard Hamiltonian is defined by, 

+V 2J rcio-"j<r' ~M /; n i ( 2 ) 

(ij)o-er' i 

where c\ a (c. ) are annihilation (creation) operators for an electron with spin a in the Wannier orbital at site i, 
1l ia — cJa-Ciff is f ne density operator for electrons of spin er, t is the hopping matrix element and \i is the chemical 
potential. The terms proportional to U and V represent respectively the on-site and nearest- neighbor interactions. 
The latter (V) vanishes for the usual Hubbard model. 

To generalize the TPSC approach, we use the same formalism as in Ref. The Green function is defined by 

= -{JVc (1)4 (2)} (4) 

where the brackets () represents a thermal average in the canonical ensemble, T T is the time-ordering operator, and 
t is the imaginary time. The space and imaginary time coordinates (ri,n) are abbreviated as (1) . The generating 
functional Z[4> a ] is given by 

Z[4> a ] = Tr [e-P H T T S((3)] (5) 

and S(f3) is defined as follows 

lnS(/3) = -yW drdT'cl(T)c ia (r')4> a (i,lr,T'). (6) 

,4.41 



The equation of motion for the Green function takes the form of the Dyson equatior 

G- 1 (l,2)=G - 1 (f,2)-^(f,2)-S ff (I,2). (7) 

The expression for EG comes from the commutator with the interaction part of the Hamiltonian so that the self-energy 
in the above equation is given by 
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S CT (1, 2) = -U (^4(1)^(1)^(1)4(3)) G-\% 2)-Vj2 (TtcIO- + a ) c °'( l + a)c CT (l)4(3), ) G~ l (3, 2) (8) 

a' ,a 

where 5 = — er and the summation on a runs over the next neighbor sites of the site 1 (at the same imaginary time). 

The first approximation is inspired by TPSC. We factorize the two-body density matrix operators as follows, 
defining numbers with overbars, such as 3, to mean a sum over the corresponding spatial indices and integration over 
the corresponding imaginary time: 

£ CT (1, 2) = UG 5 {1, 1+)G CT (1, 3)G; 1 (3, 2)^(1, 1) + V ^ G a ,(l + a, 1 + a+)G CT (l, 3)G- X (3, 2)<w(l, 1 + a) 

<j' ,a 

= US(1, 2)G ff (l, 1 + ) 5ct? (1, 1) + VS{1, 2) G °'( l + M + « + )w (1, 1 + «)• (9) 

If we set g aa '(i, j) = 1 in the above expression, one recovers the Hartree-Fock approximation for the U term, and the 
Hartree approximation for the V term. The Fock term in the near neighbor interaction V leads to a renormalization 
of the band parameters and is not as essential as in the on-site term, as discussed in our previous paper.— The 
above expression for the self-energy goes beyond the Hartree-Fock approximation because of the presence of the pair 
correlation functions g a a'{hj) defined by 

9 ™'^ 2 >- <n CT (l))(n CT ,(2)) ' (1Uj 

This quantity is related to the probability of finding one electron with spin a' on site j when another electron with 
spin er is held on site i. The presence of the pair correlation function is analog to the local field correction in the Singwi 
approach to the electron gas42 Thanks to this correction factor, we recover the exact result that the combination 
E CT (1, 3)G CT (3, 2) is equal to the potential energy when point 2 becomes 1 + (where the superscript + refers to the 
imaginary time and serves to order the corresponding field operator to the left of operators at point 1). For this 
special case then, g aa r(l, 1 + a) contains the information about the Fock contribution as well. 

We need the response functions obtained from the first functional derivative of the Green function with respect to 
the external field. Taking the functional derivative of the identity G CT (1, 3)G~/(3, 2) = 5 aa >5(l, 2) and using Dyson's 
equation Eq. ([7|), we find 



W (1,2; 3, 3) = -^^3j = GMA) tfM3>3) 0.(5,2) 



= -^G.(l,3)G g (3,2) - E^(M) g^ g#^ (M)- 



£0(6,7) 5^(3,3) 

The irreducible vertex is the first functional derivative of the self-energy respect to the Green function, as can be seen 
from the last equation. It can be obtained from our equation for the self-energy Eq. ([9]): 

Sr^rk = U6a*»5(4, 5)5(4, 6)5(5, 7). 9 ^(4, 4) + l/V 5(4, 5)5(4 + a, 6)5(5 + a, 7)<w- (4, 4 + a) 
5G CT »(6,7) ^ 

In the standard Hubbard model with V = 0, there is only one unknown functional derivative of the pair correlation 
function with respect to the Green function. Since it enters only the charge response function , it can be determined by 
making the local approximation, enforcing the Pauli principle and requiring that the fluctuation-dissipation theorem 
relating charge susceptibility to pair correlation function be satisfied^ In the present case, there are not enough sum 
rules to determine the additional functional derivatives of the pair correlation functions with respect to the Green 
function. 

In our previous paper on the extended Hubbard model^ we neglected these additional unknown functional deriva- 
tives. Here we present an approximation for these functional derivatives that improves our previous results, as can be 
judged by comparisons with QMC calculations. We would like to replace the unknown terms by functional derivatives 
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with respect to the density as fgf^^fj ~ ^f^TTp) ^(^' 4)- Unfortunately we cannot evaluate this type of functional 
derivative either, unless point 3 coincides with either 2 or 1, so we approximate the previous equation by 

5g aa > (1,2) Sg a<7 , (1,2) Sg a<J , (1,2) 

drier" (3) on CT "(l) on a »(2) 

In other words, we ignore the terms where point 3 gets further away from either 1 or 2. We will show that the last 
equation is exact at half-filling (see Eq. (|30|) ) and that it can be evaluated . We will also show in the next section by 
explicit comparisons with QMC calculations that its extension away from half-filling is a reasonable approximation. 
The last two approximations lead to the local approximation of TPSC as a special case. 

Inserting in the susceptibility Eq. (jTTJ) the irreducible vertex obtained from Eq. (|12p and our approximation for the 
functional derivatives, we find 



X™'(l;2) = -<y CTCT ,G CT (l,2)G CT (2,l) 

-E/G CT (l,3)G ff (3,l).g a5 (3,3) 



6G a (3,3+) 
*<M2,2) 



VG a (l,3)G a (3, 1) £ <W'(3, 3 + a) ^"f + g 3 + ° +) 



£/G a (l,3)G ff (3,l)n 5 (3)^ 



VG a (l,S)G a (3,l) n °'"& 



Ha- (2,2) 

8g aS {Z,Z) 6G a „(3,3+) 
Sn a „(3) 5^,(2,2) 

Sg aa ,„ (3, 3 + a) 5G a „ (3, 3+) Sg aa ,„ (3, 3 + a) SG a „ (3 + a, 3 + a+) 



a a a 



(Sn ff »(3) (S^(2,2) fe CT "(3 + a) 



,(2,2) 



(14) 



For the repulsive case, we expect that spin and charge components of the response functions will be dominant. Thus 
we need the dynamic correlation functions of the density n = n-\ + nj. and magnetization m = n-\ — nj, respectively. 
They are given by 



Xcc,m(1;2) =2[x ffff (l;2)± X aa(l;2)] 

= -2G CT (l,2)G CT (2,l) 

-2£/G CT (l,3)G ff (3,l). 9ff5 (3,3) 



4VG CT (1, 3)G CT (3, 1) ^5cc lS ,(3, 3 + a) 



ft? g (3,3+) 5G CT (3,3+) 
50s (2, 2) 50^(2,2) 

50^(3 + 0,3 + 0+) , (5G CT (3 + a,3 + a+) 



5^(2,2) 



± 



50,(2,2) 



2n[/G ff (l,3)G CT (3,l) 
4^0,(1,3)0,(3,1)^ 

a 

4nVG CT (l,3)G CT (3,l)^ 



5.9,, (3, 3) 5^(3,3) 
5n(3) 5m(3) 



5G,(3,3+) ± <SO,(3,3^ 



^(2,2) 



5ff s ,(3, 3 + a) Sg sa (3,3 + a) 
5n(3) ' 5m(3) 

5g s ,(3, 3 + a) 5g s ,(3, 3 + a) 
5n(3 + a) 5m(3 + a) 



6.(2,2) 
5G,(3,3+) ± 5G,(3,3+) 



5^(2,2) 50,(2,2) 
5G,(3 + a, 3 + a+) 5G,(3 + a, 3 + a+) 



r(2,2) 



± 



50,(2,2) 



(15) 



where the first term A in the curly brackets {A, B} is associated with the upper + sign, and the second term B is 
associated with the — sign. Also, in the above equation, 

gcc,ss — \9<T<7 + 9<ra)/2 

are the charge and spin correlation functions, while 

g scT = (sw +sw)/2 

is the symmetric combination of the spin-resolved pair correlation functions. Indeed, the definition of g cc in partially 
spin polarized system is different from the above equation and that is the reason for defining g sa . In fact g cc is 
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symmetric function under a flip of the spins of the electrons but not g SIJ . This is an important issue when we deal 
with the functional derivative of the pair correlation functions with respect to the magnetization. The functional 
derivative terms with respect to magnetization in the fourth and the last term are zero . This is due to the fact that 
<7 CT ct(1, 1) an d flw(l, 2) are even functions under exchange of a <-» a on site 1 and on site 2 respectively. The functional 
derivative with respect to density in the fifth and last term are equal (this can also be shown by following the procedure 
that we explain in the following subsection). Now we can factorize the charge and spin response functions from the 
above equation. The final form of the response functions in Fourier-Matsubara space are given by 



x°(q,^n) 



x ss (q,0 = : — x o (q ^ ))rT ; (17) 



where 



U C M = U (W, 1) + n -^f ) + & (g cc {l, 1 + a) 7 (q) + ~ ( 2 + 7(q))) , (18) 

f/ ss (q) = Ug^(l, 1) - 4V 1 + a) 7 (q) + 2n ^^f^) ■ ( 19 ) 

In this equation, 7(q) is given by 7(q) = X^ a cos (9« a )' a = l-A being the dimension of the system, w„ = 27mT 
is the Matsubara frequency and x°(q, w„) is the free response function 

o ( , f dp /°(P + f)-/°(P-f) , 9n , 

x(q,w„)=/ — : * ■ (20) 

JBZ V lUJ n - £ p+q/2 + £ p-q/2 

In the above formula v is the volume of the Brillouin zone (BZ), /°(q) = l/[l + exp((e 9 -/z)/T)] is the free momentum 
distribution function and e q = —2tJ2 a C0S (<laCt) is the free particle dispersion relation. The pair correlation functions 
are related to the instantaneous structure factors by 

Scc(ri) = 1 + - / ^[^(cO-llexpftq-r,) (21) 

ffa.(ri) = - / — [5 ss (q) - 1] cxp(zq ■ r,) (22) 

where S' CCiSS (q) = S crcr (q) ± 5^ (q) are the charge and spin component of the instantaneous structure factor, the spin 
resolved instantaneous structure factor being defined by 5 CT(T '(q) = {n a {q)n a i (q)) jn with n<j(q) the Fourier transform 
of n- la . Finally, the instantaneous structure factors are connected to the response functions (or susceptibilities) by the 
fluctuation-dissipation theorem 

S C c,ss{Oi) = ~~y^ Xcc, ss (q, Wn) (23) 



B. Particle- hole symmetry for the functional derivatives 

We are left with the evaluation of the functional derivatives. We first use the particle-hole (Lieb-Mattis) canonical 
transformation c a (i) = cxp(iQ • Ti)c'\(i) (where Q = (it, it)) in the definition of the pair correlation function. In the 
second step we will use electron-hole symmetry to evaluate the functional derivative of the pair correlation function 
at half-filling. First step: When point 1 and point 2 arc different or when a ^ a' the particle-hole transformation 
gives 

= (n CT (lK,(2)) = ([1 -<(i)][i -n'AM = l-<(l)-<,(2) + <(lK,(2) g ; g ,(l,2) 
Wl ' ' Ml)rv(2) " [1 -<(!)][! -<,(2)] 1 - <(1) - <,(2) + <(1)<,(2) 1 J 
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where g 1 is the pair correlation function of holes at the density n' . Note that n inside the angle brackets are the 
operator densities. We abbreviate average densities such as (n a (l)) by n a {l). 

The above equation may look innocuous, but it leads to a very strong constraint on the pair correlation function and 
it can be used to find the exact functional derivative of the pair correlation function at half-filling where particle-hole 
symmetry is exact. To obtain this we first take n(l) = 1 + e and m(l) = which means n'(\) = 1 — e and m'(l) = 0. 
Note that n (2) is unchanged. Using the fact that nf,J.(*) = \ n {i) i TO (*)]/ 2 one obtains 

5T<t( 1 > 2 ) = T+i ■ ( 25 ) 

4 

When there is exact electron-hole symmetry, such as at half-filling, the electron pair correlation function at density n 
is equal to the hole pair correlation function at the same density, so we deduce from the above equation that 

-P-„.M1. (26) 
The summation of last equation for two different values of a gives us one of the functional derivatives 



5n(l) 

Using the same trick we can also show that 

<%u(M) 



8n(l) 



= [l- 5cc (l,2)]. (27) 



= 2[l- 5u (l,l)]. (28) 



To evaluate the functional derivative respect to magnetization, it suffices to assume that n(l) = 1 and m(l) = e . It 
is then easy to show that 

^|^ = [1-<7 CC (1,2)]. (29) 
Using the same trick one can also show, as mentioned before, that at half filling 

£<w(i,2) 



8{n, m)(3) 



(30) 



when site 3 is different from both 1 and 2. Although there is no strict particle- hole symmetry away from half- filling, 
it is reasonable to assume that at weak to intermediate coupling the Pauli principle forces only electrons close to the 
Fermi surface to be relevant. In that case, the dispersion relation can be linearized and particle-hole symmetry can be 
safely assumed for the pair correlation functions g a(T >(i,j). The comparisons with QMC calculations presented later 
suggest that it is indeed a good approximation. 



C. Self-consistency and discussion of Extended Two-Particle Self-Consistent Approach (ETPSC) 

Substituting the above equations for the functional derivatives in the expression for the irreducible vertices U cc (q) , 
Eq.(fl8]). and U ss (q) , Eq. (fl"9)l . one can compute the charge and spin susceptibilities Eqs. ([16ll7[) given the three 
unknown pair correlation functions cw(U l)j <?cc(U 1 + a ) an d 3ss(U 1 + &)• We do not need g^ail, 1) because we 
know from the Pauli principle that it vanishes. The three unknowns can be determined self-consistently as follows. 
Use the fluctuation-dissipation theorem Eq. (f2"3"]) to obtain S cc . ss (q) from the susceptibilities and then substitute in 
the definitions Eqs. (|21l22|) of the pair correlation functions in terms of S cc . ss (c{). From this procedure, we obtain 
four equations for three unknowns since a look at the left-hand side of Eqs. (|21l22p tells us that we can compute 
5cc(0) = g cc (l,l), g ss (0) = g ss (l,l), as well as g cc (a) = g cc (l,l + a) and g ss (a) = g ss (l, 1 + a). The problem is 
thus over determined. To be consistent with the initial assumption g aa (\, 1) = from the Pauli principle, the self- 
consistent solution of these four equations have to lead to g aa (0) = g cc (0) + g ss (0) = 0. In ordinary TPSC for the 
on-site Hubbard model, one does not have a separate equation for 8g a a{\,l)/8n(\) that enters U cc . Instead it is 
determined by enforcing that g aa (0) = in the sum rules Eqs. (|21l22p . In the present extension of TPSC to near- 
neighbor interaction we have an expression for that functional derivative, Eq. (|2"8)l so in general the four equations 
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obtained from the sum-rules Eqs. (|21l22p do not satisfy the on-site Pauli principle g cc (0) + g ss (0) = g aa (0) = 0, which 
is assumed in their derivation. 

At this point there are three possible ways to proceed, (a) Keep the equation |22|) for g ss (0) and drop that for g cc (0) 
in Eq. (f2"Tj) (b) Do the reverse procedure, i.e. keep the equation (|2"Tj) for g cc (0) and drop that for g ss (0) in Eq. ({2"2")) . 
(c) In the original spirit of TPSC, do not use our estimate of the functional derivative l)/5n(l) Eq. (|28|) but 

determine the value of 

g<ra{l, 1) + n- 



6n(l) 

entering the charge vertex Eq. (fT8|) by solving the four equations for g cc (0) — ffcc(l, 1)> 5ss(0) = g S s0-> 1), ffcc(a) = 
g cc (l, 1 + a) and g ss (a) = g ss (l, 1 + a) that follow from Eqs. (|2II22|) . The latter approach, (c), reduces to TPSC when 
the nearest-neighbor interaction V vanishes and gives the best approximation in that case. With this procedure, 
9aa (0) = is still true. However, comparisons with Quantum Monte Carlo calculations also suggest that in the 
presence of V, procedure (a) is a better approximation. When V vanishes, one recovers TPSC for the spin fluctuations. 
These fluctuations are anyway dominant, while the charge fluctuations, different from TPSC, are reduced compared 
with the non- interacting case. With the latter procedure, (a), the sum rule g cc (0) + g ss (0) = that follows from the 
Pauli principle (n 2 a ) — (n a ) is not satisfied, but the equation that allows us to find U ss when V — does take into 
account the Pauli principle and gives exactly the same result as TPSC. 

Following the same reasoning as in Appendix A of Ref.|3), we also find that our approach satisfies the / sum rule 
and the Mermin- Wagner theorem. Following the lines of Appendix B of Ref.(ll), one can also develop a formula for 
the self-energy at the second step of the approximation that satisfies the relation between EG and potential energy. 
The latter is a kind of consistency relation between one- and two-particle quantities. 



III. NUMERICAL RESULTS 

A. Accuracy of the approach 

Following the procedure (a) explained at the end of the previous section, we present in the first six figures various 
results that allow us (i) to discuss some of the physics contained in the pair correlation functions and (ii) to establish 
the validity of the approach and the sources of error by comparison with benchmark Quantum Monte Carlo calculations 
and with other approximations. All the calculations are for the two-dimensional square lattice. In most cases, unless 
otherwise specified, we will consider how the physics changes when V is varied for U = 4 and n — 1 fixed. We are 
using units in which ?i, fcs , i and lattice spacing a are all taken to be unity. When the spin or charge components of 
the response functions are growing rapidly at low temperature towards very large values, we stop plotting the results. 
The last two figures of the paper contain general crossover diagrams for the various forms of charge and spin orders, 
neglecting pairing channels. 

In Figs. 1 and 2, we first discuss the case V = 0. This is necessary since our approximation (for charge fluctua- 
tions) is different from TPSC even in that case. Fig. 1 compares results of our numerical calculations with QMC, 
TPSC and FLEX. We depict the charge (right panel) and spin (left panel) components of the static (zero frequency) 
susceptibilities at n — 0.87, U = 4 and V — as a function of temperature. Consider the spin susceptibility on the 
left. As we can see, the FLEX approximation underestimates the spin susceptibility. The response function within 
TPSC stays slightly below benchmark QMC results although the TPSC result for structure functions^ (not presented 
here) are almost on top of QMC results. The results from the ETPSC method developed in this paper coincide with 
the TPSC results for the spin component. The charge susceptibility differs between TPSC and ETPSC as illustrated 
on the right panel. We do not have QMC results for the charge susceptibility but the charge structure functions for 
TPSC and QMC are quite closed although not as much as in the case of the spin structure functions. We can thus 
surmise that FLEX underestimates the charge susceptibility. In ETPSC, the charge susceptibility is larger than in 
TPSC. 

Even though the charge fluctuations are less important than the spin fluctuations when V vanishes, it is instructive 
to consider them to gauge the effect of our two approximations, namely factorization with a local field factor Eq.([9]), 
approximate and estimate of the functional derivatives in Sec. IIIBI We compare in Fig. 2 the charge component 
of the TPSC instantaneous structure factor with ETPSC, where the functional derivative terms are evaluated, and 
with our previous approach* 1 - (WOFD stands for Without Functional Derivatives) that neglects completely these 
functional derivatives. The temperatures are relatively high but they are low enough to show discrepancies between 
the various approaches (TPSC essentially coincides with QMC results in this range of physical parameters). The 
importance of the functional derivatives is obvious from all figures in the different panels since ETPSC always agrees 
better with TPSC than WOFD. Given the assumption that the pair correlation functions are functionals of the density 




FIG. 2: Comparison of the ETPSC results for the instantaneous charge structure factor with TPSC and with WOFD. 
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FIG. 3: The pair correlation functions as a function of temperature at n — 1, U = 4, for various values of the nearest-neighbor 
interaction V = 0, 0.5, 1 and 1.25. 

fEq . (fT3]) ) . our evaluation of the functional derivatives at half-filling is exact. Hence the three panels for n = 1 in Fig. 2 
are particularly instructive. They suggest that the main difference between ETPSC and TPSC is due either to the 
factorization Eq.®, or to the assumption Eq. (|13[) that the pair correlation functions are functionals of the density 
only. In other words, the difference between ETPSC and TPSC suggests the limits of the factorization assumption. 
TPSC seems to compensate the limitations of the factorization by fixing the functional derivative term with the Pauli 
principle. One can also find the degree of violation from the Pauli sum-rule by comparison of our data with TPSC. 
This is the main source of error in our calculation and it is more pronounced at larger U, lower T and near half-filling. 
The bottom-right panel shows that ETPSC is in better agreement with TPSC away from half filling, an encouraging 
result. 

To illustrate the physics associated with the nearest-neighbor interaction V, we show in Fig. 3 the variation of 
<7(tct(0) and g C c,ssip) with temperature for various values V = 0, 0.5, 1 and 1.25. We first notice that g CT ,5-(0) at 
V = has a sharp decrease around T « 0.3, which means that the probability for finding two particles at the 
same place sharply decreases with decreasing the temperature around that point. g ss (a) also becomes more negative 
around the same temperature, which means that the probability of finding of two electrons at a distance a and with 
opposite spins is higher than finding them with the same spins. These two results are a consequence of the tendency 
towards Spin-Density Wave (SDW) order at zero temperature. The decrease of g a &(0) around T w 0.3 signals then 
entry into the renormalized classical regime, as discussed in Ref. 0. In the charge sector, g cc (a) does not show 
any strong change in the low temperature limit, which means there are no strong Charge Density Wave (CDW) 
fluctuations at V — 0. We observe that this feature changes as we increase V, which shows that SDW fluctuations 
are depressed in favor of CDW fluctuations. This can easily be observed at V = 1.25, which shows that g C c(a) has a 
sharp decrease around T w 0.5 and instead ga&(0) has a sudden increase. On the other hand g ss (a) dose not show any 
decrease around that temperature. These results together show for V = 1.25 the buildup of CDW rather than SDW 
fluctuations as temperature decreases. We stress that all these low temperature behaviors are non-perturbative. The 
high temperature behavior where the quantities vary smoothly can be understood perturbatively. 

To gauge the accuracy of our approach at finite V, we now compare ETPSC results with benchmark QMC 
calculations 10 for various values of V — 0.5, 1 and 1.25. We plot in Fig. 4 the static (zero frequency) spin and 
charge response functions evaluated at Q = (it, it) as a function of temperature along with the QMC results in the 
left and right panels respectively. As is clear from the figure, there is a good agreement between our results and 
the QMC results. As V increases, the charge response is largest and for that quantity the agreement improves at 
larger V. A simple comparison of this figure with the results of our previous paper— shows the importance of the 
functional derivative terms and the accuracy of their evaluation. One also observes that the spin component of the 
response function exhibits a rapid decrease at low temperature limit when V = 1.25. An analogous behavior occurs 
at V = (not presented here) for the charge response function. This feature is a result of the Pauli principle, which 
implies that in the paramagnetic state charge and spin fluctuations cannot both increase^. The small deviation of the 
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FIG. 4: The spin (left panel) and charge (right panel) zero-frequency response functions at n = 1, U = 4, V = 0.5, 1 and 1.25. 
Points are the QMC results^ and lines are the ETPSC results. 
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FIG. 5: The spin (left panel) and charge (right panel) static structure factors at n 
are the QMC results^ and curves are the ETPSC results. 



1, U = 4, V = 0.5, 1 and 1.25. Symbols 



charge response (or instantaneous structure factor) at V = 0.5 can be removed by implementing this sum-rule^. We 
also observe that the spin component of the response functions at different values of V all tend to the same value at 
sufficiently high temperature, which simply means that g ss (l, 1 + a) and ^fo^Tj are negligible in this limit. This is 
the only region where ETPSC coincides with Random Phase Approximation (RPA) results and the only place where 
RPA is valid. The effects of the above-mentioned terms become very important at low temperature and that is one 
of the sources of failure of RPA. 

Fig. 5 shows the instantaneous spin and charge structure factors at Q = (tt,tt) and V = 0.5, 1 and 1.25 as a 
function of temperature. All the features are similar to what was mentioned in Fig. 4 except that the agreement is 
less satisfactory. This suggests that our approach with constant vertices is best for low frequency response functions. 
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FIG. 6: The density-density correlation function in term of the temperature at V = 0.5, 0.75 and 1 for (7 = 0. Symbols are 
the QMC results^ and curves are the ETPSC results. 

The structure factors contain contributions from all frequencies and we know for example that the irreducible vertices 
should become equal to the bare interaction at large frequency 4 . The ETPSC approach that we are using, described 
in procedure (a) Sec. Ill C[ gives better agreement with QMC at higher values of V. In the case of the charge structure 
factor, procedure (c) compares more favorably with QMC for small V. 

Fig. 6 shows the density-density correlation function ((n 2 ) = ((nj, + wp 2 )) as a function of temperature for V — 
0.5, 0.75 and 1 in the extreme case U = 0. The agreement with QMC is quite remarkable. In this case, we used 
procedure (b) in Sec. Ill CI that fixes gaa(0) from the charge rather than the spin component of the response functions. 
The difference with procedure (a) is small, but in this limiting case where charge fluctuations dominate, procedure 
(b) is the one that agrees best with QMC. 

B. Crossover temperatures towards renormalized classical regimes 

The Mermin- Wagner theorem prevents breaking of a continuous symmetry in two dimensions. Instead, as tem- 
perature decreases one observes a crossover towards a regime where the characteristic fluctuation frequency of a 
response function at some wave vector becomes less than temperature (in dimensionlcss units). This is the so-called 
renormalized-classical regime. In that regime, the correlation length grows exponentially until a long-range ordered 
phase is reached at T = 0. The wave vector of the renormalized classical fluctuations does not depend only on the 
non-interacting susceptibility and it does need to be guessed, it is self-determined within ETPSC. Based on previ- 
ous experience with TPSC and on numerical results of the previous section, ETPSC will fail below the crossover 
temperature when the correlation length becomes too large, but the crossover temperature itself is reliable. 

We present in Figs. 7 and 8 the crossover temperature as a function of U and V for n = 1 and n — 0.75 respectively. 
Even though we display results for negative U and V as well, other instabilities, in particular in pairing channels, can 
dominate when one of the interactions is negative. One would need to generalize ETPSC for negative interactions 
along the lines of Refs. 00, but this has not been done yet. The results for negative interactions are given for 
illustrative purposes only. 

The crossover temperatures in Figs. 7 and 8 are obtained from the condition \cc {Qx, 1y, 0)/xo(<faj Qy, 0) = 10 or 
Xssilx, %, 0)/xo(Qx, Qy, 0) = 10, where q x and q y refer to the wave vector where the response function has a maximum. 
These quantities are proportional to £ 2 , the square of the corresponding correlation length. The characteristic spin 
fluctuation frequency scales like £ -2 . When the crossover temperature is reached, £ increases so rapidly that a choice 
different from 10 does not change very much the estimate and can increase sensibly the computation time. A more 
detailed discussion on the crossover temperature can be found in Ref. ITU 

The red lines in Fig. 7 separate three different fluctuation regimes, namely (it, it) SDW (on the left), (it, ir) CDW (on 
the right) and phase separation (PS) (this is defined by the growth of the zero frequency charge response function at 
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FIG. 7: The crossover temperature at n — 1. The red lines separate (tt, 7r) SDW (left), (n,n) CDW (right) and PS (on top). 
Results for negative values of the interaction strength are given for illustrative purposes only. Pairing instabilities may dominate 
in these cases. 

q x = and q y = 0). It is noteworthy that the boundary between the SDW and CDW crossover is closely approximated 
by U w 4V (for positive U and V), which is in good agreement with previous results^ i 10 ' 29 i 36 ' 37 i 38 . Theory that includes 
collective modes beyond mean-field also shows deviations from U = 4Y at weak to intermediate coupling^ Since 
there are four neighbors which, through V, favor charge order while U favors spin order, the crossover at U « 4V is 
expected. In fact, the result U — zV, where z is the number of nearest-neighbors, seems to hold quite generally in 
various dimensions at T — 0,l i 19 i 39 Slightly away from half-filling (n — 0.9) and with a small second-neighbor hopping, 
the MF approximation 35 predicts a coexistence region for SDW, CDW and ferromagnetism in the large U and V region 
of parameter space. This is likely to be an artefact due to an inaccurate estimation of the spin fluctuations in MF 
when V ^ 0. One can see the weakness of this approach by setting g ss — and ignoring its functional derivative in 
the equation for the spin susceptibility Eq. (|17[) . These terms are the main source of the suppression of the SDW so 
that one obviously overestimates the spin effects by ignoring these terms when V is large. One can easily notice that 
the charge and spin fluctuations appear at relatively large temperature over most of parameter space. In addition, 
negative U greatly increases charge fluctuations. These will compete with pairing fluctuations that we have not taken 
into account here. At positive U and V, d— wave superconductivity is not likely to occur in regions other than that 
which is nearly black in the figure. 

Fig. 8 shows the crossover diagram for n = 0.75. The fluctuations outside the red lines are commensurate. The 
commensurate CDW persists away from half-filling 36 while the SDW fluctuations are greatly reduced in the U > 
0, V > regime. The fluctuations are incommensurate in the narrow area between the green and red lines. For U and 
V inside the green curve, the system is paramagnetic, which means that neither the spin nor charge response functions 
can satisfy the above mentioned criteria even at very low temperature (T = 0.01). Finally the blue lines separate the 
region of parameter space where the spin fluctuations are larger than charge fluctuations (left bottom) from the region 
where the reverse is true (right top). Commensurate fluctuations have the highest crossover temperature. The next 
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FIG. 8: The crossover temperature for n = 0.75. Results for negative values of the interaction strength are given for illustrative 
purposes only. Pairing instabilities may dominate in these cases. In the region outside the red line, commensurate (1) SDW 
(n, tt), CDW (-7T, 71") and PS (0, 0) fluctuations dominate. (The regions of interest are similar to those observed in Fig 7 at n = 1, 
except that SDW become much less important). In the area surrounded by the red and green lines, the renormalized classical 
fluctuations are incommensurate ISDW in (2) and ICDW in (3). The area inside the green line boundary is paramagnetic. 
The blue line separates the regions where the spin (left and bottom) and charge (right and top) components of the response 
functions are dominant. 
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most important fluctuations are incommensurate. They exist next to commensurate fluctuations but their crossover 
temperatures are rather small. 

IV. CONCLUSION AND SUMMARY 

In this paper, we extended the non-perturbative TPSC approach to the case where we add to the standard Hubbard 
model a nearest- neighbor interaction V. We verified the validity of this Extended Two-Particle Self-Consistent approach 
(ETPSC) by comparisons with Quantum Monte Carlo calculations. ETPSC applies up to intermediate coupling and 
it satisfies conservation laws and the Mermin- Wagner theorem. It also includes renormalization of the irreducible 
vertices by quantum fluctuations (Kanamori-Bruckner) . 

In the original TPSC, the evaluation of the irreducible vertices requires one pair-correlation function and its func- 
tional derivative. These two quantities are evaluated by imposing the Pauli principle on two sum rules that lead 
to self-consistency conditions. In the presence of V, a similar procedure for the irreducible vertices allows us to 
determine the needed pair correlation functions self-consistentlyii, but not the functional derivatives of these pair 
correlation functions. One approach is to neglect the functional derivatives, in the spirit of the Singwi approach^ to 
the electron gas, as we did earlier—. The new idea of the present paper is to use particle-hole symmetry to evaluate 
the missing functional derivatives. When V is finite, we found that using the resulting estimate of the functional 
derivatives and dropping one of the sum rules between spin and charge fluctuations implied by the Pauli principle 
gives better agreement with benchmark QMC than that obtained by neglecting the functional derivatives^. Although 
particle-hole symmetry strictly applies only to the nearest-neighbor hopping model at half-filling, one expects that 
the approximation will not be too bad away from half-filling. The results for n = 0.8 on the bottom right panel of 
Fig. 2 are encouraging in this respect. 

For illustrative purposes, we also presented the crossover diagram of the two-dimensional extended Hubbard model 
on the square lattice for n = 1 and n — 0.75 and a range of values of V and U. Incommensurate fluctuations appear 
away from half-filling, but usually they do so at temperatures much lower than commensurate fluctuations. We 
highlighted the differences with mean-field calculations^. The crossover temperature indicates when the fluctuations 
associated with a self-consistently determined wave vector become large, or more specifically when these fluctuations 
enter the renormalized classical regime. Based on earlier arguments^, we expect that a pseudogap appears in the 
renormalized classical regime whenever the wave vector of the fluctuations can scatter electrons from one side to the 
other side of the Fermi surface. 

The methodology developed in the present paper can be easily applied to different lattices and higher dimension. 
It can be extended also to longer-range interactions. Note also that the functional derivatives introduce just two 
independent unknowns, namely ^^jy^ and ^f^fe^ = S9 sm(i)^ ' ^ there were independent ways of estimating the 
compressibility and the uniform spin susceptibility, that may give, with the Pauli principle, new sum-rules that could 
be used to estimate these functional derivatives. That question and that of multiband models are for future research. 
We also plan to apply ETPSC to specific compounds, such as the cobaltates, where charge fluctuations are important. 
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